home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / RUNGKUT.LBR / REACT.FOR < prev    next >
Text File  |  1986-04-11  |  2KB  |  79 lines
  1. $NOFLOATCALLS
  2. $NODEBUG
  3. $STORAGE:2
  4. c************************************************************************
  5.  
  6.     implicit double precision (a-h,p-z)
  7.     dimension y(2),work(34),icom(4)
  8.     external freact
  9.     common /reacts/ifeval,Da,delta,beta,Hw,Tw
  10.  
  11.     open(2,file= ' ',status= 'new')
  12.     ifeval=0
  13.     icom(1)=0
  14.     icom(2)=0
  15.     icom(3)=0
  16.     write(*,*)  'Wall Temp.=, Reactant inlet Temp=, htc='
  17.     read(*,*) Tw,Tr,U
  18.     Tw=Tw/Tr
  19.     U=U/1000.0
  20.     write(*,*) ' imeth=, tola=, tolr='
  21.     read(*,*) imeth,tola,tolr
  22.  
  23. c****   evaluate constants in the equations
  24.     Da=2.d0*5.d0/3.d0
  25.     beta=0.03d0*1.d04/1.2d0/1.d0/Tr
  26.     delta=1.d3/8.3144d0/Tr
  27.     Hw=2.d0*U*2.d0/0.1d0/1.2d0/1.d0/3.d0
  28. c****
  29.     hstart=0.01d0
  30.     neqn=2
  31.     x0=0.d0
  32.     xb=0.d0
  33.     y(1)=1.d0
  34.     y(2)=1.d0
  35.     conc=y(1)*0.03
  36.     temp=y(2)*Tr
  37.     write(2,30)xa,y(1),y(2),hstart
  38.     do 20 j=1,10
  39.        xa=xb
  40.        xb=0.1*dble(j)+x0
  41.        call runkut(xa,y,xb,neqn,tola,tolr,hstart,work,
  42.      &                 imeth,ierror,icom,freact)
  43.  
  44.        conc=y(1)*0.03
  45.        temp=y(2)*Tr
  46.        if(ierror.GT.1)then
  47.          write(2,30)xb,y(1),y(2),hstart
  48.          write(2,*)' ERROR-Problem too stiff or is discontinous'
  49.          close(2)
  50.          stop
  51.        else
  52.          write(2,30)xb,y(1),y(2),hstart
  53.        end if
  54. 20      continue
  55.     if(icom(4).GT.0) write(2,*) ' Severe round-off error possible'
  56.     write(2,*) ' Number of function evaluations = ',ifeval
  57.     close (2)
  58.     stop
  59. 30      format(1x,f11.2,1x,d15.6,1x,f15.4,1x,d15.6)
  60.     end
  61. c***********************************************************************
  62.  
  63. c       user supplied subroutine containing the system of first
  64. c       order ordinary initial value differential equations
  65.  
  66.     subroutine freact(x,y,yprime,neqn)
  67.  
  68.     implicit double precision (a-h,p-z)
  69.     dimension y(neqn),yprime(neqn)
  70.     common /reacts/ifeval,Da,delta,beta,Hw,Tw
  71.  
  72.     yprime(1)= -Da*y(1)*dexp(delta*(1.d0-1.d0/y(2)))
  73.     yprime(2)= beta*Da*y(1)*dexp(delta*(1.d0-1.d0/y(2)))
  74.      &             -Hw*(y(2)-Tw)
  75.     ifeval=ifeval+1
  76.  
  77.     return
  78.     end
  79.